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Ultracold atomic gases are, to a very good approximation, described by pairwise zero-range inter¬ 
actions. This paper demonstrates that A^-body systems with two-body zero-range interactions can 
be treated reliably and efficiently by the finite temperature and ground state path integral Monte 
Carlo approaches, using the exact two-body propagator for zero-range interactions in the pair prod- 
net approximation. Harmonically trapped one- and three-dimensional systems are considered. A 
new propagator for the harmonically trapped two-body system with infinitely strong zero-range 
interaction, which may also have applications in real time evolution schemes, is presented. 


I. INTRODUCTION 

Systems with two-body zero-range interactions consti¬ 
tute important models in physics. Although realistic two- 
body interactions do typically have a finite range, results 
for systems with zero-range interactions provide a start¬ 
ing point for understanding complicated few- and many- 
body dynamics. In 1934 [1], Fermi used the zero-range 
model in quantum mechanical calculations to explain the 
scattering of slow neutrons off bound hydrogen atoms. 
Nowadays, the two-body contact interaction is discussed 
in elementary quantum texts [2]. It has, e.g., been used 
to gain insights into the correlations of molecules, such as 
and H 2 , and to model atom-laser interactions [3-5]. 

In the 60s [6-10], zero-range interactions were used ex¬ 
tensively to model strongly-interacting one-dimensional 
systems at zero and finite temperature. Many of these 
models are relevant to electronic systems where the 
screening of the long-range Coulomb interactions leads 
to effectively short-range interactions [11]. More recently, 
ultracold atomic gases interacting through two-body van 
der Waals potentials have been, in the low temperature 
regime, modeled successfully using zero-range interac¬ 
tions [12-14]. One-, two-, and three-dimensional systems 
have been considered. 

While zero-range interactions have been at the heart 
of a great number of discoveries, including the Efimov ef¬ 
fect [15-17], their incorporation into numerical schemes 
is not always straightforward. Loosely speaking, the 
challenge in using zero-range interactions in numerical 
schemes that work with continuous spatial coordinates 
stems from the fact that we, in general, do not know 
how to incorporate the boundary conditions implied by 
the zero-range potential into numerical approaches at the 
four- and higher-body level. 

This paper discusses an approach that allows for the 
use of zero-range potentials in many-body simulations. 
We work in position space and consider a system with 
fixed number of particles. We develop a scheme to in¬ 
corporate pairwise zero-range interactions into e~'^^ di¬ 
rectly, where H is the system Hamiltonian. The quan¬ 
tity is of fundamental importance. If r is iden¬ 

tified with l/iksT), where ks and T denote the Boltz¬ 


mann constant and temperature, respectively, then 
is the density matrix for the system at finite tempera¬ 
ture. Knowing the density matrix, the thermodynamic 
properties can be calculated. If, on the other hand, r is 
identified with it/h, where t denotes the real time, then 
6“’’^ can be interpreted as the real time propagator and 
be used to calculate dynamic properties. Throughout 
this paper, we refer to r as imaginary time, keeping in 
mind that r carries units of 1/energy and that r can be 
associated with inverse temperature or real time. 

The remainder of this paper is organized as follows. 
Section II reviews the pair product approximation, which 
relates the many-body propagator to the two-body prop¬ 
agator. Section III derives the two-body propagator for 
various systems with zero-range interactions. Sections IV 
and V demonstrate that the two-body zero-range prop¬ 
agators yield reliable results if used in one- and three- 
dimensional path integral Monte Carlo (PIMC) [18-20] 
and path integral ground state (PIGS) [18, 21-24] sim¬ 
ulations of trapped iV-atom systems. The performance 
and implementation details will be discussed. While the 
free-space zero-range propagators have been reported in 
the literature [25-28] , the zero-range propagators for the 
harmonically trapped system with infinite coupling con¬ 
stant are, to the best of our knowledge, new. Finally, 
Sec. VI concludes. 


II. V-BODY DENSITY MATRIX 

We consider N particles with mass mj and position 
vector Tj (j = 1,... ,7V) interacting via a sum of zero- 
range potentials with interaction strength g. The Hamil¬ 
tonian H of the system can be written as 

N N 

H = + ( 1 ) 

j=i j<k 

where denotes the non-interacting single-particle 

Hamiltonian of the jth particle and Vjk the two-body 
potential between the jth and kth particle. In the fol¬ 
lowing, it will be convenient to separate the Hamilto¬ 
nian Hjk, where Hjk = -|- Vjk, of atoms j 
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and k into relative and center of mass pieces, Hjk = 
Hjf + , where depends on the relative vec¬ 
tor rjk, and on the center of mass vector 

J'ifc = I’i “ and -|- mkYk)/{mj -|- m^). 

Below, the non-interacting two-particle system will serve 
as a reference system and we define H^f. = 
and 77]^ =77™'-°+ 7?;,-. 

The A^-particle density matrix ptot(R-, in posi¬ 

tion space can be written as 


the importance of 7V-body {N > 2) correlations. If r is 
identified with l/[kBT), the pair product approximation 
is limited to high temperature. In this case, the pair 
product approximation is analogous to a virial expansion 
that includes the second-order but not the third-order 
virial coefficient [30]. 

III. TWO-BODY RELATIVE DENSITY MATRIX 


Ptot(R,R';r) = (R|e-""|R'), (2) 

where R = (ri,...,rAr) and R' = (r'^,...,r^) collec¬ 
tively denote the coordinates of the 7V-particle system. 
For sufficiently small r, ptot(R, R^; t) can be constructed 
using the pair-product approximation [18], 

/ N 

Ptot(R,R';T) ~ 


In the following, we consider one- and three- 
dimensional systems, without and with external har¬ 
monic confinement, and discuss the evaluation of the 
relative density matrix for zero-range interactions. For 
notational simplicity, we leave off the subscripts j and k 
throughout this section, i.e., we denote the relative dis¬ 
tance vector by r for the three-dimensional system and 
X for the one-dimensional system, respectively, and the 
relative part of the two-body Hamiltonian by 77''®*. 



I , (3) 

j<k ) 


where p’'®*(rjfc, r'^; r) denotes the normalized pair density 
matrix. 


p''®*(rjfc,r'fc;T) = 


fl''®*(rj-fc,r'-fc;T) 


(4) 


and p''®*(rjfe, r'^; r) and p''®*’°(rjfc, r',.; r) the relative den¬ 
sity matrices of the interacting and non-interacting two- 
body systems. 


p-'(r,,,r' ;t) 






(5) 


A. One-dimensional system 

The complete set of bound and continuum states of 
Tf'®* is spanned by ipn with eigen energies En and ipk 
with energies where fi denotes the reduced 

two-body mass and k the relative scattering wave vector. 
If the ipn and ij^k are normalized according to 

J(x) da: = 6nn' (8) 

and 

Jil^l{x)'ipk'ix)dx = S{k - k'), (9) 


and 


/®*’°(rjfc,r'fe;T) = (rjk 


^rel,0 

e“ J" 




( 6 ) 


In Eq. (3), p®P(rj,r';T) denotes the single-particle den¬ 
sity matrix. 


'Ir ■ r'■ 






(7) 


The key idea behind Eq. (3) is that the one- and two-body 
density matrices can, often times, be calculated analyt¬ 
ically. Indeed, the non-interacting propagator is known 
in the literature both for the free-space and harmonically 
trapped systems [18, 29]. Moreover, the eigen energies 
and eigen states of the Hamiltonian 77'|* have, for a class 
of two-body interactions, compact expressions, which en¬ 
ables the analytical evaluation of the relative two-body 
density matrix in certain cases (see Sec. HI). 

It should be noted that the pair product approximation 
is only valid in the small r limit since it does not account 
for three- and higher-body correlations. For the real time 
dynamics, this means that the time step is limited by 


then the relative density matrix p'^®*(a:, x': t) can be writ¬ 
ten as [29] 

/"'(a:, x'; r) = ^ tjj*^ix)e~'^^’-''ipn{x') + 

n 

poo 

/ ( 10 ) 

Jo 

Free-space system: The relative Hamiltonian for the 
free-space system with zero-range interaction can be writ¬ 
ten as 

where g denotes the coupling constant of the ^-function 
potential. For positive g, the Hamiltonian given in 
Eq. (11) does not support a bound state and the cor¬ 
responding energy spectrum is continuous. The sym¬ 
metric and anti-symmetric scattering states with energy 
h?k'^/{2p) read [31] 

= -^sm{k\x\-^S{k)) 

VTT 


( 12 ) 
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and 


i^kix) = ^sm{kx), (13) 

y/TT 

respectively; S{k) = arctan [fi^fc/(5/i)] is the phase shift. 
For negative g, the Hamiltonian additionally supports a 
bound state with symmetric wave function 

^o(^) = (14) 

and energy —g^/i/(2/i^). Integrating over the symmet¬ 
ric and anti-symmetric scattering states, and adding, for 
negative g, the additional bound state, one finds the nor¬ 
malized relative density matrix [25-28] , 


Trapped system: For two particles in a harmonic trap, 
the system Hamiltonian reads 


= - 


2pL dx"^ 


-gd{x) 


9 2 
- UUJ X , 


(18) 


where w denotes the angular trapping frequency. The 
energy spectrum of is discrete and the eigen ener¬ 
gies and eigen functions are known analytically in com¬ 
pact form [35]. These solutions can be used to evaluate 
Eq. (10) numerically. The corresponding relative non¬ 
interacting density matrix reads 


rel,0 


1 - 1/2 


FiD,trap(2^>2:';T) = [27r sinh(rfcj)ah, 

(x^ -I- x'^) cosh(T/kt;) — 2xx' 


exp 


2 sinh(T/Ej)a^ 


ho 


(19) 


-rel / / N 1 f At(a;x' -I- lxx'l)\ 

PiD,free(a;,a: ;r) = 1 - exp (- — -1 


Ti/iT g 
2 h 


erfc(u) exp(u^), (15) 


where u = g, (Jx'] -I- [x] -|- gr) /yj2 and erfc is the 
complementary error function. We emphasize that 
Eq. (15) holds for positive and negative g. The corre¬ 
sponding relative non-interacting density matrix reads 



g{x-x'f 
2 tK^ 


(16) 

The free-space propagator given in Eq. (15) was em¬ 
ployed in a PIMC study of the harmonically trapped 
spin-polarized two-component Fermi gas with negative 
9 [32]. 

For large \g\, u approaches ^gxI2gjh and, using 
limu_).oo ^/^^^^(u) exp(M^) = 1, Eq. (15) reduces to 


PiD.free(a;)a:';T) 


1 — exp for xx' > 0 

0 for xx' < 0. 

(17) 


Equation (17) suggests that the relative coordinate does 
not change sign during the imaginary time evolution. 
Since the interaction strength is infinitely strong, the 
two particles fully reflect during any scattering process, 
i.e., the transmission coefficient is zero. This means that 
the initial particle ordering remains unchanged during 
the time evolution. This is a direct consequence of the 
Bose-Fermi duality of one-dimensional systems [6, 33, 34]. 
Specifically, the phase shift of the symmetric wave func¬ 
tion given in Eq. (12) goes to zero when jgj —>■ oo, imply¬ 
ing that the symmetric wave functions coincide, except 
for an overall sgn(x) factor, with the anti-symmetric scat¬ 
tering wave functions of non-interacting fermions. The 
implications of the Bose-Fermi duality for Monte Carlo 
simulations of A'^-body systems with infinite g is discussed 
in Sec. IV. 


where Oho denotes the harmonic oscillator length, Oho = 
yjhjijiw). For fixed r and finite 5, one can then tabu¬ 
late Pid\rap(2^i 2;'; r) for discrete x and x' using Eq. (10) 
and use a two-dimensional interpolation during the /V- 
body simulation. The infinite sum in Eq. (10) can be 
truncated by omitting terms with n > nmaxj where rimax 
is chosen such that the Boltzmann factor fulfills the in¬ 
equality <C The value of rimax depends on 

the time step: smaller r require larger rimax- 

For infinite g^ we were able to derive a compact analyt¬ 
ical expression for PiD_trap(®i *^5 5 to infinity, 

the probability distribution of each even state coincides 
with that of an odd state, i.e., the system is fermionized. 
The complete set of even and odd eigen states for g = 00 
can be written as 


and 


fl.{x) = (/>n(la;l) 

(20) 

''t^ni.x) = 4>n{.x), 

(21) 


where (j)n{x) is the non-interacting harmonic oscillator 
wave function. 


= (^2”n!aho) "'"/(^“S°)iF„(x/aho), (22) 


Hn{x) denotes the Hermite polynomial of order n, and 
n takes the values 1, 3, 5, 7,.... The corresponding ener¬ 
gies are En = {n + l/2)huj for both the symmetric and 
anti-symmetric states, i.e., each energy level is two-fold 
degenerate. Using Eqs. (20) and (21) in Eq. (10) and 
evaluating the infinite sum analytically, we find 


—rel / f \ 
PlD,trap(^i^ 



for xx' > 0 

for xx' < 0. 
(23) 


For rhuj 1, i.e., when the trap energy scale is much 
smaller than 1/r, the trap propagator [Eq. (23)] equals 
the free-space propagator [Eq. (17)]. 

To test the one-dimensional propagators for infinite g, 
we consider the Hamiltonian given in Eq. (18) and pre¬ 
pare an initial state using a linearly discretized spatial 
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grid. Our aim is to determine the ground state wave func¬ 
tion and energy by imaginary time propagation. Two 
approaches are used. First, the initial state is propa¬ 
gated using the exact trap propagators [see Eqs. (23) and 
(19)]. In this case, the error originates solely from the 
discretization of the spatial degree of freedom; indeed, 
we find that the energy approaches the exact ground 
state energy quadratically with decreasing grid spacing 
5x. Second, the initial state is propagated using the 
free-space propagator [see Eqs. (17) and (16)]. We ap¬ 
ply the Trotter formula [36] and move half of the trap 
potential to the left and half to the right of the free- 
space Hamiltonian. This is known as the primitive ap¬ 
proximation [18], which is expected to yield a quadratic 
time step error since the trap potential does not com¬ 
mute with the free-space Hamiltonian. The error is found 
to scale quadratically with both the time step and the 
grid spacing. For r = (50fiw)“^, 5x = v^Oho/dO, and 

^max = ^min = where Xniin ^ ^ ^ ^max; 

we obtain energies that deviate by 2.4 x and 

2 X for the free-space propagator and the trap 

propagator, respectively, from the exact ground state en¬ 
ergy of 3/ia;/2. 


B. Three-dimensional system 

Because the s-wave zero-range potential is spherically 
symmetric, the relative orbital angular momentum op¬ 
erator commutes with the relative Hamiltonian. Corre¬ 
spondingly, we label the bound states ipnim with eigen 
energies E^i and the continuum states il^kim with ener¬ 
gies / (2/r) by the relative orbital angular momentum 
quantum number I and the projection quantum number 
m. If the ipnim and tpkim are normalized according to 

J' (^d) 

and 

J'^klm (r) dr = (5(fc - k')Sw 5mm’ , (25) 

the relative density matrix ^^'(r, r'; r) can be written 
as [29] 

nlm 

_ nOO 

^ / rkimir)e-^’^'^'/^^^Hkim{v')dk. (26) 

Im -^0 

Free-space system: The Hamiltonian of the three- 
dimensional system in free space reads 

= + (27) 

p or 

where ag is the s-wave scattering length. The second 
term on the right hand side of Eq. (27) is the regularized 


two-body zero-range pseudopotential [37]. The I > 0 
continuum states read 


i>kim{r) = -kji{kr)Yim{r), (28) 

V TT 

where the ji and Yim denote spherical Bessel functions of 
the first kind and spherical harmonics, respectively. The 
continuum states for the s-wave channel read 

V’fcoo(i’) = — sin(A:r -f Sg{k)), (29) 

y/ZTTr 

where Ss{k) = arctan(—a^fc) is the s-wave phase shift. 
Eor positive Ug, there exists an s-wave bound state with 
eigen function 'i/'ooo(r) = l/-\/27raI^exp(—r/og) and 
eigen energy —h^/{2pa1). As is evident from the above 
eigen states, only the s-wave states are affected by the 
interactions. Thus, we construct the relative interacting 
density matrix by writing the non-interacting rel¬ 
ative density matrix and subtracting from it the 

non-interacting s-wave contribution and adding to it the 
s-wave contribution for finite Og. 

Eor negative Og, there exist only continuum states and 
the density matrix can be expressed as 

PrD.free(l', r'; r) =P 3 D(Le(r, r'; t) + 

Th^k^ 1 

e 2 a. ^—-[sin(fcr-I-5g(/c)) sin(fcr'-b (5g(fc)) — 
ZTT'^rr 

sin(fcr) sin(fcr')] dfc, (30) 



where the non-interacting relative density matrix reads 

orDi..(r.rV) = (31) 


The integral in Eq. (30) can be done analytically [38] and 
the normalized relative density matrix reads [27, 28] 


-rpl A / X 1 f 

PsD.freelr, r ; t) = I -b - - exp - 

prr' \ 


prr'{I -b cos 6*) 
1 t?t 


I -I- \ ®^fc(u) exp(?;^) 

Qj Q y ZtjJj 


(32) 


where cosO = r • r'/{rr') and v = [r-br' — 
rh^/ (pog)]/ \j2Tk?jp. Adding the bound state contri¬ 
bution to Eq. (30) [see the first term on the right hand 
side of Eq. (26)] for positive Og, one finds Eq. (32), i.e., 
the same propagator as for negative Og [28]. 

Eor [ogl = oo, Eq. (32) simplifies to 




P3D,free(r, y'-t) = I + — exp - 


prr' 


prr'(l -b COS 0 ) 

n^T 


(33) 

This propagator was recently used in a proof-of-principle 
diffusion Monte Carlo study of the homogeneous two- 
component Fermi gas at unitarity with zero-range inter¬ 
actions [39]. 
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Trapped system: The Hamiltonian for two particles in 
a spherically symmetric harmonic trap with s-wave scat¬ 
tering length Us reads 


_rpl ^ r-r2 ^ 2 2 ft fig cfSl / \ ^ /o 

2/i 2 /i or 


The non-interacting relative density matrix reads 

PrD.trapCr^i-V) = Oho^ pTTsinh(rSa;)]x 
(r2 _|_ j./2^ cosh(rSw) — 2r • r' 


exp 


2sinh(T/ia;)a^Q 


(35) 


Similar to the free-space case, the relative interacting 
density matrix is obtained from the non-interacting den¬ 
sity matrix with the difference of the s-wave eigen states 
and energies of the interacting and non-interacting sys¬ 
tems added. For finite Og, we were not able to evaluate 
the inhnite sum analytically. Because of the rotational 
invariance, the infinite sum depends only on r and r' (and 
not the direction of the vectors r and r'), allowing for an 
efficient tabulation of the reduced relative density matrix. 
For infinitely large Og, we hnd an analytical expression. 
In this case, the bound state wave functions that are af¬ 
fected by t he (5- function interaction can be written as 
■\/2(()n('c)/V47rr2, where the are defined in Eq. (22) 
with X replaced by r and n = 0,2,4,... The relative 
two-body density matrix reads 


rel / / \ rel.O / / \ , 

P3D.trap(i’, r ; t) = P3D(t,,,p(r, r ; r) -b 

E (36) 

n—0 


where the even n terms in the sum over n come from the 
s-wave states that are affected by the interactions and 
the odd n terms come from the s-wave states of the non¬ 
interacting system. Performing the infinite sum, we find 
for the normalized relative density matrix 


P3D.trap(l',r';T) = 1-b 

^2 

sinh(T?ia;) exp 


rr'(l + cosO) 
sinh(rfiw) 


(37) 


Setting rhu = 0, Eq. (37) reduces to Eq. (33), i.e., to the 
corresponding free-space expression. 

Equations (17) and (33) show that the one- and three- 
dimensional free-space propagators for systems with in¬ 
finitely large i5-function strength are characterized by the 
length y'rWJJi, which is proportional to the de Broglie 
wave length. The trap propagators for infinitely large 
coupling constant [see Eqs. (23) and (37)], in contrast, 
are characterized additionally by the harmonic oscilla¬ 
tor length. For finite interaction strength, the coupling 
constant defines a second length scale for the free-space 
system and a third length scale for the trapped system. 



N 


Figure 1. (Color online) PIMC results for N harmonically 
trapped distinguishable one-dimensional particles with two- 
body zero-range interactions of infinite strength at tempera¬ 
ture T — fujj/kB- (a) The symbols show the energy obtained 
by the PIMC approach as a function of the number of par¬ 
ticles N. For comparison, the dotted line shows the exact 
thermally averaged energy, (b) Symbols show the energy dif¬ 
ference AS between the PIMC energies and the exact results. 
As a reference, the dotted line shows the AS = 0 curve. In 
(a) and (b), the circles and squares are calculated using 8 and 
128 time slices, respectively. 


IV. ONE DIMENSIONAL TESTS 

This section incorporates the trapped two-body prop¬ 
agator into PIMC calculations for one-dimensional N- 
particle systems with pairwise zero-range interactions. 
We find that the conventional PIMC sampling ap¬ 
proaches [18, 19] yield an efficient and robust description 
of one-dimensional systems with two-body zero-range in¬ 
teractions. 

As a first example, we consider N distinguishable har¬ 
monically trapped particles with mass m in one spatial 
dimension with pairwise zero-range interactions of infi¬ 
nite strength. The V-particle system with infinitely large 
interaction strength is unique in that the particle statis¬ 
tics becomes irrelevant for local observables. For exam¬ 
ple, the energy is the same for N identical bosons, N 
identical fermions and N distinguishable particles at any 
temperature provided all particles interact via two-body 
zero-range interactions. We employ the zero-range trap 
propagator together with the single-particle trap propa¬ 
gator. Symbols in Fig. 1(a) show the PIMC energy for 
N distinguishable particles at temperature ksT = huj. 
Circles and squares are for simulations with imaginary 
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time step tHuj = 1/8 and 1/128, respectively. For com¬ 
parison, the dotted line is calculated directly from the 
partition function of N non-interacting fermions. Figure 
1(b) shows the energy difference AE between the sim¬ 
ulation and the analytical results. It can be seen that 
the calculations for the larger time step (circles in Fig. 1) 
exhibit a systematic time step error, which is found to 
scale quadratically with the time step r for fixed N and 
to originate from the pair product approximation. Since 
we include the two-body correlations exactly, the leading 
order error is expected to come from three-body corre¬ 
lations. Indeed, for relatively small fixed r (t“^ around 
and varying N, we find that the error AE scales 
approximately linearly with the number of triples, sug¬ 
gesting that the error is dominated by three-body corre¬ 
lations with sub-leading contributions arising from four- 
body correlations. As the number of particles N or the 
time step r increase, four- and higher-body correlations 
become more important. This error analysis suggests 
that an improved propagator could be obtained if the 
three-body problem could be solved analytically. We 
note that the performance of the zero-range trap propa¬ 
gator for the system with infinite g is quite similar to that 
of the free-space propagator using the second- or fourth- 
order Trotter decomposition. The reason is that the error 
is dominated by three- and higher-body correlations. 

We now discuss that the simulations need to be mod¬ 
ified to treat N identical bosons or fermions with pair¬ 
wise zero-range interactions of infinite strength. Equa¬ 
tions (17) and (23) indicate that the paths for any two 
particles cannot cross. This implies that the permute 
move, implemented following the approach discussed in 
Ref. [23], yields a zero acceptance probability. This is 
consistent with our finite g simulations for N identical 
bosons. As we change g for otherwise hxed simulation 
parameters from small positive to large positive values, 
the probability of sampling the identity permutation ap¬ 
proaches 1. The fact that particle permutations are al¬ 
ways rejected, causes two problems for the infinite g simu¬ 
lations. First, since the particle ordering does not change, 
the single particle density for the first particle differs from 
that of the second particle, and so on. To calculate the 
single particle density p{x) of, e.g., the N identical boson 
system, one can average the single particle density pj (x) 

of the jth particle over all j, p{x) = Al~^Y^^=iPj{^)- 
An analogous average can be performed for other lo¬ 
cal (closed paths) structural properties. Second, the 
simulation of open paths, which allow for the calcula¬ 
tion of off-diagonal long-range order, requires that the 
sampling scheme be modihed since open paths do al¬ 
low for permutations. The two-particle density matrix 
/9({xi,X 2 },{x 2 ,Xj};r), e.g., is finite if xi < X 2 < x'^. 
The calculation of non-local observables is beyond the 
scope of the present paper. 

As a next example, we apply the zero-range trap prop¬ 
agator to iV = 2 and 3 identical bosons in a harmonic 
trap with g = fi^/(-\/2/raho)- For the PIMC calculations, 
we tabulate the density matrix for the time step of in¬ 


terest and interpolate during the simulation. For small 
number of particles, we expect the zero-range trap prop¬ 
agator to work well even for a large time step and we 
use rhuj = 1/2. For ksT = hLijl32, we obtain an energy 
of 1.3067(l)fiw and 2.3880(l)/ia; for = 2 and 3, re¬ 
spectively. The temperature is so low that the system is 
essentially in the ground state. For comparison, we deter¬ 
mined the ground state energy using the transcendental 
equation from Ref. [35] and by solving the Lippmann- 
Schwinger equation [40]. The resulting ground state en¬ 
ergies [1.306746fcu and 2.3880(l)fiw for iV = 2 and 3, 
respectively] agree within error bars with the PIMC re¬ 
sults. 

To demonstrate that the PIMC simulations describe 
the short-distance behavior of systems with zero-range 
interactions correctly, we analyze the pair distribution 
function Pi 2 {x), which is normalized to f_^Pi 2 (x) dx = 
1, for N = 2 and 3 identical bosons with finite g. To 
start with, we derive the short-distance properties of the 
pair distribution function Pi 2 for N identical bosons with 
zero-range interactions. Using the Hellmann-Feynman 
theorem, the partial derivative of the energy with respect 
to g can be related to the pair distribution function at 
X = 0 [41], 


Fi2(0) 


2 dE 
N{N-1) dg ■ 


(38) 


Note that Eq. (38) is the one-dimensional analog of 
equating the three-dimensional adiabatic and pair rela¬ 
tions [42]. Second, from the Bethe-Peierls boundary con¬ 
dition of the fV-body wave function 4* (the derivatives 
are taken while all other coordinates are kept hxed), 


94- 


dd! 



dxjk 

Xjfc-i>0+ 

dxjk 

Xjff; ^0 

A2 


(39) 


one can derive that the slope of the pair distribution func¬ 
tion at any temperature satishes 


^X^12(x) ,^^ 0 + 

- 


-40 ■ 

Eor identical 

bosons. 

9Pi2(x)/9x 

(40) 

and 

dPi2{x)/dx\^^^_ 

have the 

same magnitude 

but 

opposite signs. 





The dashed and dotted lines in Eig. 2 show Pi 2 (x) 
obtained from our PIMC simulation for = 2 and 
3, respectively. Eor comparison, the solid lines are ob¬ 
tained using Eqs. (38) and (40). The values of dE/dg 
are obtained through independent energy calculations us¬ 
ing the techniques of Refs. [35, 40]. We hnd Pi2(0) = 
0.3266002/aho and 0.308245(2)/aho for A^ = 2 and 3, re¬ 
spectively. Our PIMC results agree well with the solid 
lines in the small |x| regime, demonstrating that the 
PIMC approach describes the short-range behavior ac¬ 
curately. 
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Figure 2. (Color online) PIMC results for harmonically 
trapped one-dimensional bosons interacting through two- 
body zero-range interactions with coupling constant g — 
fi^/(%/2/raho) at temperature T = /ki;/(32fcs). The dashed 
and dotted lines show the pair distribution function obtained 
by the PIMC approach for N = 2 and 3 bosons, respectively. 
For comparison, the solid lines show the asymptotic short- 
range behavior obtained by alternative means (see text). 



Figure 3. (Color online) Scaled pair distribution functions 
47rPi2(r)r^ for two distinguishable particles of mass m in 
a harmonic trap at kBT/{hijj) — 1. The solid and dashed 
lines are for two particles with infinitely large s-wave scatter¬ 
ing length interacting through the zero-range potential and 
a Gaussian potential with effective range re « 0.0861aho, re¬ 
spectively. For comparison, the dotted line is for the non¬ 
interacting system. 


V. THREE DIMENSIONAL TESTS 

The pair distribution function of the non-interacting 
three-dimensional system is, unlike that of the non¬ 
interacting one-dimensional system, zero at vanishing in¬ 
terparticle distance. This fact leads, as we discuss now, to 
non-ergodic behavior unless the traditional path integral 
sampling methods are complemented by an additional 
move. To motivate the introduction of this new “pair 
distance move”, we consider the two-particle system. 

Solid and dashed lines in Fig. 3 show the scaled pair 
distribution functions for two distinguishable particles 
with infinitely large s-wave scattering length in a three- 
dimensional harmonic trap at ksT/{fko) = I interacting 


through a zero-range potential and a finite-range Gaus¬ 
sian potential with effective range Ri O.OSGlohoj re¬ 
spectively. The pair distribution function Pi 2 {r) is nor¬ 
malized according to dvr Pi 2 {r)r'^ dr = 1. The pair 
distribution functions for the finite-range and zero-range 
potentials nearly coincide for large r, but differ for small 
r. The pair distribution function for the Gaussian po¬ 
tential drops to 0 as r approaches 0 while that for the 
zero-range potential approaches a non-zero value. 

PIMG and PIGS simulations typically use the first 
term on the right hand side of Eq. (3) as the “prior distri¬ 
bution” and the second term as the “correction”. This is 
suitable for iV-body systems with two-body finite-range 
potentials for which the scaled pair distribution func¬ 
tion is, as that of the non-interacting system (see the 
dotted line in Fig. 3 for a two-body example), zero at 
r = 0. Since the prior distribution has zero probability 
at r = 0 , the pair distribution function of the system 
with zero-range interactions is not properly sampled if 
standard sampling schemes are used. Ergodicity is vio¬ 
lated at r = 0 and the probability to sample the region 
near r = 0 is small. Moreover, the correction term [see 
Eqs. (37) and (33)] diverges as r or r' go to 0. This means 
that the probability to sample configurations with r r; 0 
is small. However, if such a configuration is chosen, there 
is a small probability to accept a new configuration with 
much larger r, i.e., the correlation length is large. Sim¬ 
ilarly, if one uses the naive uniform distribution for the 
prior distribution, i.e., if one proposes a move for which 
all Cartesian coordinates differ by 5x from the current 
configuration, where 5x is a random value between — Aa; 
and Ax, the problems discussed above still exist. 

To remedy the problems that arise if the standard sam¬ 
pling approaches are used, we introduce a “pair distance 
move” for which the prior distribution scales as 1 /r^ in 
the pair distances. First, particles j and k and a bead I 
are chosen (the coordinates for the 1 th bead are collec¬ 
tively denoted by Rj) and the distance rjk = \rjk\ and 
the direction Yjk are calculated. A new vector rjk^new 
that lies along Tjk or —is proposed, r^fc^new = 

The quantity e is written as e = rjk -|- Sr, where Sr is 
obtained by choosing uniformly from —Ar to Ar. If the 
weight w, 

. r, Ptot (R/ — 1 ; R/,neW; ^)ptot (R/-t-l j R/,new ; , 

w = mm 1 ,-—-—---—-—- -—5 - , 

Ptot (R/-1, Ri; Tjptot (R;-|-1, Ri; T)r A 

(41) 

is larger than a uniform random number between 0 and 
1, then the proposed move is accepted. Otherwise, it 
is rejected and the old configuration is kept. The value 
of Ar is adjusted such that about 50% of the proposed 
moves are accepted. It can be easily proven that detailed 
balance is fulfilled. Our PIMC calculations for the two- 
body system with zero-range interactions show that the 
“pair distance move” significantly improves the sampling. 
Without this move, the short-range behavior of the pair 
distribution function has a long correlation length, which 
increases with decreasing r. With this move, the small 





r behavior is described accurately. The move described 
here is related to the compression-dilation move intro¬ 
duced in Ref. [43]. Few details were given in Ref. [43] 
and no comparison with that approach is made in this 
paper. 

We now demonstrate that the outlined sampling 
scheme provides a reliable description of Bose systems at 
unitarity, which have attracted a great deal of attention 
recently experimentally and theoretically [43-50]. While 
the properties of unitary Fermi systems with zero-range 
interactions are fully determined by the s-wave scattering 
length [12, 14, 51-54] those of Bose systems additionally 
depend on a three-body parameter [15, 17]. Specifically, 
if the two-body interactions are modeled by zero-range 
potentials, then a three-body regulator is needed to pre¬ 
vent the Thomas collapse of the 7V-boson {N > 3) sys¬ 
tem [17, 55]. Here, we utilize a purely repulsive three- 
body potential of the form 

VMm) = (42) 

^jki 

where Rjki denotes the three-body hyperradius, Rjki = 
-I- Tji + rki)/3. In the A^-boson system, each of the 
iV(7V — 1)(7V — 2)/6 triples feels the regulator, i.e., the 
term J2j<k<i ^shiRjki) is added to the Hamiltonian with 
pairwise zero-range interactions. In the absence of an 
external trap, the zero temperature three-body ground 
state energy iftrimer of the unitary system is set by the 
Cq coefficient. The corresponding length scale is 1 /k, 
where k = m\Etriiner \/h is the binding momentum. 

Our goal is to determine the ground state properties 
of self-bound A^-boson droplets at unitarity in the ab¬ 
sence of an external confinement. In the context of the 
present paper, it would seem that our goal could be read¬ 
ily achieved using the PIGS approach. It turns out, how¬ 
ever, that without a good initial trial wave function, the 
number of time slices needed to converge the calculations 
is rather large, making the simulations computationally 
quite expensive. Instead, one might consider performing 
PIMC calculations at various temperatures and extrapo¬ 
lating to the zero temperature limit. This approach also 
turns out to be computationally extensive. Our simula¬ 
tions pursue an alternative approach, in which the scat¬ 
tering states of the system are discretized in such a way 
that the relative ground state energy iSduster of the N- 
body cluster is much larger than the energy scale intro¬ 
duced by the discretization. We utilize a spherically sym¬ 
metric harmonic trap and adjust the trapping frequency 
such that [Flciusterl fkuj. Simulations are then per¬ 
formed at a temperature where the Bose droplet is in the 
ground-state dominated liquid-phase [43, 56], where the 
finite temperature introduces center-of-mass excitations 
but not excitations of the relative degrees of freedom. 
The temperature Ttr at which excitations of the relative 
degrees of freedom become relevant can be estimated us¬ 
ing the “combined model” introduced in Ref. [56] . As we 
show now, this approach allows for a fairly robust deter¬ 
mination of the Ai-boson properties at zero temperature. 



Figure 4. (Color online) Scaled pair distribution function 
47rPi2(r)r^ for three identical harmonically trapped three- 
dimensional bosons with two-body zero-range interactions 
with infinitely large s-wave scattering length and repulsive 
l/P® three-body potential. The solid line and squares are 
calculated by the zero-temperature PIGS approach and the 
PIMC approach at T — ttw/ks- For comparison, the cir¬ 
cles show the scaled pair distribution function obtained by 
sampling the exact ground state density using the Metropolis 
algorithm. 

We set the trap energy huj to 0.27|£'trimer| (Tiitrimer is 
the ground state energy of the three-boson system in free 
space) and the temperature to tvjj/ks- These parameters 
provide a good compromise: First, the temperature is 
sufficiently low that finite temperature effects are negli¬ 
gible (i.e., T < Ttr for the N considered below, N = 3—9) 
and high enough that convergence can be reached with 
the computational resources available to us. Second, the 
size of the A^-boson system is smaller than the harmonic 
oscillator length such that structural properties such as 
the pair distribution function are largely unaffected by 
the external confinement for A^ > 5. 

Our path integral simulations use the two-body zero- 
range trap propagator. The repulsive three-body poten¬ 
tials are treated using the Trotter formula. In the second- 
order scheme, half of the sum of the three-body potentials 
is moved to the left and half to the right of the Hamilto¬ 
nian H that accounts for the two-body interactions and 
the external confinement. In the fourth-order scheme, a 
more involved decomposition is used [57, 58] . In addition 
to the standard moves and the “pair distance move”, we 
implement a move that updates the center-of-mass coor¬ 
dinates. The introduction of this center-of-mass move is 
motivated by the fact that the relative degrees of freedom 
are expected to be, to a good approximation, “frozen” in 
the ground state while low-energy center-of-mass excita¬ 
tions are allowed. Indeed, the center-of-mass energy is 
given by E'c.m. = SHui coth{huj/{2kBT))/2, which evalu¬ 
ates to 3.24593/ia; for T = fiuj/kB, indicating that center- 
of-mass excitations cannot be neglected. 

The squares in Fig. 4 show the pair distribution func¬ 
tion calculated by the PIMC approach for three identical 
bosons at T = huj/kB- For comparison, the solid line 
and the circles show zero-temperature results. The solid 
line is calculated using the PIGS approach with a trial 
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Figure 5. (Color online) PIMC energies for N = 5 har¬ 
monically trapped three-dimensional bosons with two-body 
zero-range interactions with infinite scattering length and 
and repulsive 1/77® three-body interaction at temperature 
T = hujjkB as a function of the time step r. The cir¬ 
cles (lower-lying data set) and squares (higher-lying data set) 
show the energy obtained using the second- and fourth-order 
scheme, respectively. The error bands are obtained by fitting 
the data for different r intervals. 


wave function that coincides with the exact ground state 
wave function [59] while the circles are calculated by sam¬ 
pling the exact zero-temperature ground state density 
using the Metropolis algorithm. The agreement between 
the three sets of calculations is very good, demonstrat¬ 
ing (i) that excitations of the relative degrees of freedom 
are negligible at the temperature considered and (ii) that 
the path integral approaches accurately resolve the short- 
range behavior of the pair distribution function. The pair 
distribution functions shown in Fig. 4 are affected by the 
trap, i.e., they move to larger r as the trap frequency uj is 
reduced. The reason is that huj is only about four times 
smaller than |i?trimer| • The magnitude of the 7V-boson en¬ 
ergy i?ciuster increases rapidly with N [60-62], implying 
that the trap effects decrease quickly with increasing N^ 
thus allowing us to extract the free-space energy iSciuster 
from the finite-temperature trap energies i?sim- 


Figure 6. (Color online) Free-space A'^-boson ground state 
energy ^cluster as a function of N for infinitely large two- 
body s-wave scattering length. The circles with errorbars are 
extracted from our PIMC simulations. The dotted line shows 
the energies reported in Ref. [60]. 


second- and fourth-order scheme, respectively. The free- 
space energy i?ciuster is then obtained by subtracting the 
center-of-mass energy, Flciuster = - E^.m.- 

The squares in Fig. 6 show Flciuster for = 5 — 9. 
The corresponding energies E^i^ are obtained using the 
fourth-order scheme with tHoj Ri 0.000122. As can be 
seen from Fig. 5, this energy lies within the extrapolated 
T = 0 errorbands for = 5. Figure 6 scales the energy 
^cluster by the corresponding zero-temperature free-space 
trimer energy Atrimer calculated by the hyperspherical 
coordinate approach. For comparison, the dotted line 
shows the Af-boson ground state energies from Ref. [60] 
for a finite-range two-body Gaussian potential with in¬ 
finitely large s-wave scattering length and a hardcore 
three-body regulator. While the model interactions dif¬ 
fer, the agreement between the two sets of calculations 
is good, providing further support for the (approximate) 
universality of Af-boson droplets. More detailed compar¬ 
isons will be presented elsewhere [63]. 


Symbols in Fig. 5 exemplarily show our PIMC en¬ 
ergies Esim for the five-boson system at T = hw/kB 
as a function of the time step r. Circles and squares 
are obtained using the second- and fourth-order schemes 
(see earlier discussion), respectively, to treat the term 
exp(—^ TV 3 b(Rjfc;))- The statistical errors are 
smaller than the symbol size. The fourth-order results 
display, as expected, a smaller time step dependence than 
the second-order results and are well described by a func¬ 
tion of the form cq -b C 2 r^ -I- C 4 r'^, whereas the second- 
order results are described by a function of the form 
Co -b C 2 T^. The presence of the term for the fourth- 
order results is due to the fact that the pair product 
approximation neglects three- and higher-body correla¬ 
tions (see also Sec. IV). The shaded regions in Fig. 5 
show errorbands obtained by fitting the two sets of PIMC 
energies for different r intervals. The errorbars of the ex¬ 
trapolated zero time step energies are found to overlap. 
We find Asim = —37.0(1.2)fia; and —36.2(I.0)Iia; for the 


VI. CONCLUSION 

This paper discussed how to treat zero-range two-body 
interactions in V-body Monte Carlo simulations. We 
showed that the incorporation of the exact two-body 
zero-range propagator via the pair product approxima¬ 
tion allows for an accurate description of paradigmatic 
strongly-interacting one- and three-dimensional model 
Hamiltonian. 

An important aspect of the studies presented is that 
the strength of the contact interaction requires no renor¬ 
malization since the simulations are performed using 
the regularized two-body zero-range pseudopotential and 
continuous spatial coordinates in an unrestricted Hilbert 
space. The fact that the interaction strength does not 
need to be renormalized distinguishes the simulations 
presented in this paper from lattice approaches [64-66] 
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and from configuration interaction approaches [67-69]. 

The developments presented in this paper open a num¬ 
ber of possibilities. The use of two-body zero-range in¬ 
teractions, e.g., provides direct access to the two-body 
Tan contact [42] , without extrapolation to the zero-range 
limit. The two-body Tan contact is defined for systems 
with two-body zero-range interactions. It relates dis¬ 
tinct physical observables such as the large momentum 
tail and aspects of the radio frequency spectrum, and 
has attracted a great deal of theoretical [41, 70-75] and 
experimental [46, 76-79] interest. From the computa¬ 
tional perspective, the adiabatic relation, which involves 
the change of the energy with the scattering length, and 
the pair relation [42, 80], which gives the probability 
of finding two particles at the same position, are most 


convenient. Earlier work applied these relations to sys¬ 
tems with finite-range interactions and extrapolated to 
the zero-range limit. Using our zero-range propagators, 
these relations can be used directly for the determination 
of the Tan contact, eliminating the extrapolation step. 
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